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ABSTRACT 
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There  is  great  potential  for  the  use  of  satellite-  or  aircraft-home  synthetic  aperture 
radar  (SAR)  imagery  for  analysis  of  meteorological  systems.  However,  the  interpretation 
of  SAR  imagery  requires  increased  understanding  of  processes  that  change  the  state  of  the 
sea  surface.  Wind  stress  that  results  from  kilometer-scale  boundary  layer  spanning  eddies 
(ELSE)  produces  sea-surface  stress  patterns  responsible  for  some  of  this  sea-surface 
variability.  In  this  study,  a  Galerkin  model  of  these  ELSE,  which  take  the  form  of  two- 
dimensional  rolls  or  three-dimensional  convective  cells,  is  modified  with  new  boundary 
conditions  to  study  their  effect  on  sea-surface  stress. 

Previous  investigations  using  spectral  models  have  for  simplicity  been  formulated 
using  lower  boundary  conditions  in  which  the  stress  and  heat  flux  vanish.  However  these 
boundary  conditions  do  not  allow  the  study  of  stress  variability  at  the  sea  surface.  The 
marine  atmospheric  boundary  layer  model  used  here  allows  subgrid-scale  heat  and 
momentum  fluxes  at  the  lower  boundary  and  has  been  shown  in  a  previous  study  to 
capture  much  of  the  qualitative  behavior  of  ELSE  circulations.  However,  the  model 
solution  in  this  study  had  a  maximum  vertical  velocity  at  the  lower  boundary  that 
prevented  the  fluxes  from  being  physically  correct,  and  the  solution  did  not  equilibrate 
owing  to  an  unknown  energy  source.  The  goal  of  this  thesis  is  to  investigate  hypotheses 
concerning  how  these  problems  can  be  solved.  These  hypotheses  are  derived  from  an 
energetics  analysis  of  the  model  equations  that  leads  to  the  introduction  of  new  lower 
boundary  conditions  to  control  the  boundary  energy  terms.  The  model  is  then  applied  in  a 
case  study  using  Hi-Res  2  data  as  initial  conditions  to  determine  the  velocity  fields  at  the 
lower  boundary  that  produce  kilometer-scale  sea-surface  stress  variability  and  to  compare 
this  pattern  with  that  found  on  an  ERS-1  SAR  image. 
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The  case  study  results  demonstrate  that  the  model  with  the  new  boundary 
conditions  does  indeed  capture  the  spatial  organization  of  BSLE  and  correctly  reproduces 
an  expected  three-dimensional  circulation.  A  linear  analysis  using  a  small-shear,  light- 
wind  profile  and  a  -2  °C  air-sea  temperature  difference  produces  preferred  aspect  ratios 

whose  magnitudes  agree  well  with  those  derived  from  the  SAR  image.  When  these 
parameter  values  are  input  into  the  model,  a  quickly  equilibrating  solution,  with  total 
energy  that  is  fairly  constant,  is  obtained  as  expected  from  earlier  studies.  Also  as 
expected,  the  wind  profile  and  temperature  profile  modification  coefficients  become 
steady,  while  the  ELSE  circulation  and  temperature  coefficients  have  a  period  matching 
that  given  by  the  linear  analysis.  Cross-sectional  contours  of  the  streamfunctions  in  the 
across-wind  and  along-wind  directions  show  that  the  circulation  fills  the  domain,  with 
wavenumber  one  dominant,  but  with  evidence  that  other  wavenumbers  contribute  as  well. 
The  cross-sectional  results  also  show  that  the  temperature  perturbation  has  a  dominant 
wavenumber  two  pattern,  with  a  maximum  at  the  lower  boundary  where  the  thermal 
forcing  is  introduced.  These  streamfunction  and  temperature  contours  have  patterns  very 
similar  to  those  given  by  much  higher-resolution  large-eddy  simulation. 

The  most  important  result  of  the  case  study  concerns  the  BLSE-induced  sea- 
surface  stress  pattern.  At  the  lower  boundary,  this  pattern  is  obtained  from  the  vector  sum 
of  the  model-calculated  winds  and  the  background  mean  wind,  all  expressed  relative  to 
the  ocean  current.  Via  the  standard  drag  law  parameterization,  the  total  horizontal  wind 
speed  is  proportional  to  the  square-root  of  the  sea-surface  stress  magnitude.  Quite 
pleasingly,  the  stress  pattern  on  this  planview  matches  the  orientation  and  spacing  of  the 
stress  pattern  found  on  the  SAR  image.  This  result  shows  that  the  model  is  now  ready  to 
be  used  in  further  studies  of  how  ELSE  affect  the  sea-surface  stress  variability  in  different 
boundary  layer  configurations. 
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There  is  great  potential  for  the  use  of  satellite-  and  aircraft-borne  synthetic 
aperture  radar  (SAR)  as  meteorological  and  oceanographic  observation  systems  (Vesecky 
and  Stewart  1982).  The  variability  of  stress  on  the  sea  surface  produces  distinct  patterns 
in  SAR  imagery.  Secondary  circulations  within  the  marine  atmospheric  boundary  layer 
(MABL)  can  modulate  the  sea-surface  wave  field  and  so  can  produce  these  discernible 
signatures.  Evidence  for  these  secondary  circulations  has  been  found  by  Thompson  et  al. 
(1983),  Gerling  (1986),  and  Alpers  and  Briimmer  (1994),  who  studied  the  SAR  signatures 
produced  by  roll  circulations  in  the  MABL.  By  correctly  interpreting  the  signatures  in  the 
imagery,  we  can  infer  the  atmospheric  conditions  in  the  ocean  surface  layer  that  lead  to 
alternating  regions  of  converging  and  diverging  air  at  the  surface  (Gerling  1986;  Alpers 
and  Briimmer  1994).  Thompson  et  al.  (1983)  noted  that  we  must  understand  these  effects 
theoretically  before  employing  SAR  as  remote  sensors  of  ocean  surface  winds.  We 
model  these  circulations  and  use  the  model  to  study  their  effects  on  the  sea-surface  stress 
patterns  through  atmospheric  modeling  of  the  MABL  flow. 

Quasi-linear  and  cellular  patterns,  with  spacing  on  the  order  of  a  kilometer,  are 
commonly  observed  in  SAR  imagery  (e.g.,  Visecky  and  Stewart  1982).  There  are  two 
types  of  kilometer-scale  circulations  that  affect  sea-surface  stress  —  roll  circulations  and 
convective  cells.  The  first  type  is  a  quasi-two-dimensional  feature,  while  the  second  type 
is  a  fully  three-dimensional  flow  that  may  form  independently  of  the  rolls  in  less  windy 
conditions  (Woodcock  1940,1975;  Deardorff,  1976).  Etling  and  Brown  (1993)  have 
reviewed  three  fundamentally  different  mechanisms  for  the  formation  of  rolls  and  cells  in 
the  boundary  layer.  The  first  two,  the  inflection  point  and  parallel  instabilities,  are  related 
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to  the  background  wind  shear.  For  the  first  mechanism  there  must  be  an  inflection  point 
in  the  wind  profile  and  for  the  second,  the  time  scale  must  be  long  enough  that  the 
Coriolis  force  becomes  important.  The  third  mechanism  is  thermal  instability  in  which 
the  driving  force  is  an  air/sea  temperature  difference.  This  thermal  mode  produces 
convection  very  much  like  Rayleigh-Benard  convection  observed  in  the  laboratory. 
Following  Sikora  and  Young  (1993),  we  call  these  rolls  and  cells  boundary  layer 
spanning  eddies  (ELSE).  The  ELSE  vertical  scale  is  the  boundary  layer  depth,  and  they 
transfer  momentum  and  heat  down  to  the  top  of  the  surface  layer.  From  there,  small  scale 
eddies  are  responsible  for  the  final  transfer  of  momentum  and  heat  to  the  sea  surface. 

A  number  of  people  have  studied  ELSE  and  have  modeled  them  using  nonlinear 
dynamical  systems  (e.g.,  Shirer  1986;  Laufersweiler  and  Shirer  1989,  1995;  Haack  and 
Shirer  1992).  These  studies  have  sought  to  identify  the  dominant  forcing  mechanisms 
and  characteristic  spatial  and  temporal  responses  of  these  boundary  layer  circulations  as 
the  forcing  rate  varies.  These  nonlinear,  Galerkin  models,  which  are  based  on  the 
Eoussinesq  system  of  equations  and  use  sinusoidal  basis  functions,  provide  wind  profile 
modifications,  preferred  orientation  angles  and  wavelengths,  and  economy  in  the  study  of 
bifurcation  and  stability  properties.  These  models  also  yield  profiles  of  the  vertical  fluxes 
of  horizontal  momentum  and  heat.  For  simplicity,  the  lower  boundary  conditions  used  in 
these  models  are  rigid  and  stress-free.  Since  we  are  studying  the  effects  of  these  ELSE  on 
stressing  the  sea  surface,  these  boundary  conditions  are  not  appropriate  for  this  study.  We 
need  a  model  that  allows  the  heat  and  momentum  fluxes  pass  through  the  surface  layer 
and  penetrate  the  lower  boundary. 

Zuccarello  (1994)  has  developed  a  nonlinear,  Eoussinesq,  three-dimensional 
boundary  layer  model  that  allows  nonzero  subgrid-scale  heat  and  momentum  fluxes  at  the 
lower  boundary.  He  uses  similarity  theory  to  specify  the  boundary  conditions.  The 
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bottom  of  the  domain  is  set  at  10  m  above  the  sea  surface,  where  the  lower  boundary 
conditions  are  chosen  to  represent  wind  and  temperature  values  typically  observed  at  this 
height  in  the  surface  layer.  Application  of  the  stability-modified  similarity  relationships 
leads  to  constant  forcing  parameters  that  are  used  in  the  specification  of  the  lower 
boundary  conditions.  These  boundary  conditions  allow  for  a  more  realistic  incorporation 
of  thermal  forcing,  which  in  earlier  models  (e.g.,  Haack  and  Shirer,  1992)  is  spread 
throughout  the  entire  vertical  domain  and  represented  by  the  Rayleigh  number  Ra.  In 
reality,  the  thermal  forcing  at  the  lower  boundary  should  drive  the  system  and  the  value  of 
the  traditional  Rayleigh  number  should  be  inferred  from  the  model  solutions.  Finally,  and 
for  simplicity,  the  upper  boundary  conditions  are  rigid,  stress-free,  and  perfectly 
conducting  as  in  the  earlier  models. 

The  vertical  basis  functions  for  the  model  expansions  are  found  by  solving  a  more 
realistic  eigenfunction  problem  based  on  the  new  boundary  conditions  (Zuccarello,  1994). 
These  basis  functions  are  expressed  in  terms  of  exponential  and  trigonometric  functions. 
The  horizontal  basis  functions  are  composed  of  trigonometric  functions  and  are  cyclic. 

The  spectral  model,  derived  from  the  Boussinesq  system  using  these  expansions,  is 
incorporated  into  the  FORTRAN  program  that  does  the  numerical  integration  of  the 
model  and  subsequent  simulation  of  the  ELSE  circulations. 

Zuccarello  (1994)  got  fairly  good  results  in  a  simple  case  study.  The  case  study 
results  show  the  model  properly  captures  the  spatial  organization  of  the  roll  circulations. 
The  results  also  revealed  that  the  temporal  periodicity  of  the  system  is  consistent  with 
both  physical  and  numerical  analyses  of  the  solution  based  purely  on  advection  by  the 
cross-roll  mean  wind.  The  system  of  equations  and  many  of  the  physical  assumptions 
made  are  appropriate  for  modeling  ELSE  in  a  marine  environment.  However,  Zuccarello 
discovered  several  problems  during  his  case  study.  The  first  problem  is  an  unknown 
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energy  source  that  prevents  the  model  from  equilibrating  during  integration.  The  second 
problem  is  with  the  vertical  velocity.  In  the  updrafts,  the  maximum  vertical  velocity  is  at 
the  lower  boundary.  This  forces  the  maximum  momentum  fluxes  to  be  in  the  lowest  200 
m  of  the  boundary  layer,  instead  of  in  the  middle  of  the  boundary  as  typically  observed. 

The  focus  of  this  study  is  to  continue  Zuccarello’s  (1994)  work  in  studying  the 
effects  of  the  ELSE  on  sea  surface  stress  variability.  We  begin  with  an  energetic  analysis 
of  the  original  system  of  equations  and  the  model  equations  to  locate  and  identify  how  to 
control  the  energy  sources  and  sinks.  The  analysis  reveals  a  large  number  of  energy 
boundary  terms  in  the  model.  The  addition  of  some  new  boundary  conditions  eliminates 
most  of  these  boundary  terms  and  controls  the  vertical  velocity  at  the  lower  boundary.  In 
Chapter  2  we  review  the  model,  describe  the  energetic  analysis,  and  discuss  changes  we 
make  in  the  model.  In  Chapter  3,  we  look  at  some  case  studies  based  on  observational 
data  from  the  second  High  Resolution  Remote  Sensing  Experiment  (Hi-Res).  In  Chapter 
4,  we  conclude  by  reviewing  the  improvements  in  the  model  and  the  case  study  results, 
and  we  look  at  future  areas  to  study  with  this  model. 


Chapter  2 

MODEL  DESCRIPTION  AND  MODIFICATIONS 
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We  use  a  spectral  model  developed  by  Zuccarello  (1994)  to  study  boundary  layer 
spanning  eddies  (ELSE).  Zuccarello  finds  the  model  accurately  captures  the  spatial 
organization  of  roll  circulations  and  the  temporal  periodicity  of  the  system,  as  given  by 
advection  by  the  constant  cross-roll  mean  wind  component.  However,  the  model  also  has 
some  problems,  and  Zuccarello  suggests  that  an  energetics  analysis  be  performed  to  find 
the  unknown  energy  source  preventing  the  model  from  equilibrating  during  integration. 
He  also  recommends  an  examination  of  the  boundary  conditions  to  determine  why  the 
vertical  velocity  is  too  large  at  the  lower  boundary. 

We  extend  the  work  done  by  Zuccarello  by  examining  the  model  development  and 
performing  a  diagnostic  energy  analysis  to  identify  all  the  energy  sources  and  sinks  in  the 
system.  Our  results  indicate  a  number  of  boundary  terms  as  possible,  unrealistic  energy 
sources  or  sinks.  However,  the  analysis  also  indicates  that  introduction  of  additional 
boundary  conditions  will  eliminate  all  of  these  unrealistic  boundary  terms.  We  then 
modify  the  model  to  include  these  boundary  conditions  and  so  better  control  the  system 
energy. 

2.1.  Model  Description 

The  model  description  that  follows  is  a  summary  of  the  model  development  chapter 
of  Zuccarello  (1994).  For  a  complete  description  of  the  model  development,  we  direct 
the  reader  to  Zuccarello  (1994). 
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The  convective  nature  of  boundary  layer  roll  circulations  allows  the  use  of  a  model 
based  on  the  shallow  Boussinesq  system.  This  model  represents  roll  circulations  by 
perturbations  superimposed  on  a  basic  state  that  is  time-independent,  isopycnic, 
hydrostatic,  and  horizontally  moving  (Shirer  1986).  We  assume  that  the  perturbations 
possess  a  three-dimensional  structure,  and  so  we  consider  both  mean-wind-perpendicular 

and  mean-wind-parallel  variations.  The  horizontal  domain  is  infinite  and  cyclically 
continuous,  and  bounded  by  the  planes  jc  =  0, y  =  0  and  x=  L^,y  =  L^.,  where  and  L^, 

are  the  domain  wavelength  components.  If  the  domain  contained  a  single  diagonal  roll, 
then  it  would  have  wavelength  L={\l +  1/T^.)  .  The  vertical  domain  ranges  from 

z  =  hi^g  to  z  =  Zj-,  where  /i^g  is  the  height  of  the  lower  boundary  and  Zj-  is  the  cloud  top 
height  or  inversion  base  (Fig.  2.1).  Only  the  part  of  the  circulation  above  is 
specifically  modeled;  the  part  below  is  parameterized  via  boundary  conditions  based 

surface  layer  similarity  theory.  In  dimensionless  form  the  domain  is  defined  by 
0<x*<l,  0</<l,and0<z‘<l. 


Fig.  2.1.  A  schematic  cross-section  of  the  marine  atmospheric  boundary  layer,  in  which  the  surface  layer  is 
greatly  exaggerated.  In  the  surface  layer  the  temperature  decreases  and  the  wind  increases  with  height. 


The  basic  state,  known  as  the  conductive  state  and  denoted  by  subscript  q,  is  given 
by 


Vo(z)  =  V(z)  =  t/(z)i+y(z)j 

(2.1) 

Po  ~  Poo 

(2.2) 

T,(z)=T„-A^T{z/z,) 

(2.3) 

Ptt(z)  =  pi„-PMSZ. 

(2.4) 

where  subscript  qo  denotes  the  value  at  the  sea  surface,  U[z)  and  V[z)  are  the 
components  of  the  mean  horizontal  wind,  and  A^Tjzr  is  a  constant  lapse  rate  given  by 
the  temperature  difference  A^T  between  the  inversion  base  and  the  sea  surface.  The 
three-dimensional  secondary  circulation,  denoted  by  primes,  is  given  by 


u'{x,y,z,t)  =  u{x,y,z,t)-U{z) 

(2.5) 

v'{x,y,z,t)  =  v{x,y,z,t)-V{z) 

(2.6) 

w'{x,y,z,t)  =  w{x,y,z,t) 

(2.7) 

p'{x,y,z,t)  =  p{x,y,z,t)-p,. 

(2.8) 

T'{x,  y,z,t)  =  T{x,y,z,t)-T^  (z) 

(2.9) 

p\x,y,z,t)  =  p{x,y,z,t)-  pXz)- 

(2.10) 

In  the  Boussinesq  system,  compression  is  considered  negligible.  Other  approximations 
include  writing  the  equation  of  state  in  a  normalized  form  and  linearizing  the  pressure 
gradient  force  (Shirer  1987). 
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2.1.1.  The  Partial  Differential  Equations 

The  Boussinesq  system  of  equations  for  the  perturbation  temperature  T'  and  the 
perturbation  velocity  vector  v'  takes  the  form  (e.g.,  Shirer  1986) 


+ V  •  vr + v'  •  vr +KV'r 

dt  zj 

(2.11) 

+  V«Vv'  +  w'^^  +  v'*Vv'=  ^  Vp'  fkxy'  +  g 

T 

— k+vV'v' 

(2.12) 

T„ 

V*v'  =  0, 

(2.13) 

where  k  is  the  constant  eddy  thermometric  conductivity  and  v  is  the  constant  eddy 
viscosity.  A  dimensionless  form  of  the  equations  is  used  based  on  the  definitions 


,  V  K 

= - 

“  ^LB 

(2.14) 

,  W  K 

w  = - 

(2.15) 

X 

11 

X 

♦ 

(2.16) 

II 

(2.17) 

z  —  z  {zj  )  +  h^g 

(2.18) 

^  i^T  '"^Lb) 

K 

(2.19) 

rvKT;. 

g{Zj 

(2.20) 

P'= 

[zt  h^g)L^ 

(2.21) 

vW=|vfe)|v{2-) 

(2.22) 

p  * 

^  /K 
(^r  “  ^LB  ) 
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(2.23) 


These  dimensionless  equations  take  the  form 


ci  T* 

—  =  V^J*  -  Re\*  •  Vr  +  Raw*  -  v*  •  VT*  (2.24) 

ot 

1^  =  +  PV^v*  -  /*Pk  X  V*  +  P  J*k 

^  *  (2.25) 

-(PeV*  +  V*)*  Vv*  -Pew*^^, 

Vxv*=0  (2.26) 


where 


f7  ^  ^  ^  1 

ax  dy  dz 


(2.27) 


The  dimensionless  forms  lead  to  two  forcing  parameters  in  (2.24)-(2.25).  The 
Reynolds  number  Re  is  given  by 

Re  =  |V(zr)|fe-/.„)/K  (2.28) 

and  represents  dynamic  forcing  by  the  mean  wind  shear  of  the  basic  state.  The  Rayleigh 
number  Pa  is  given  by 


_  S^zT{zj 

vkToo 


(2.29) 
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and  represents  thermodynamic  forcing  across  the  domain. 

Three  other  dimensionless  parameters  appear  in  the  system.  The  eddy  Prandtl 
number  is  defined  as 


P=v/k, 


The  roll  aspect  ratios,  and  a,, ,  are  defined  as 


a,.  = 


(^r  ^LB ) 

4 

{Zr-hLs) 


(2.30) 


(2.31) 

(2.32) 


The  vector  V*(z*)  represents  the  dimensionless  mean  wind  profile  for  which  |v*(l)|  =  1 . 

Zuccarello  transforms  the  momentum  equation  (2.25)  by  taking  the  curl  to  form  a 
vorticity  equation,  thereby  eliminating  pressure.  The  resulting  equation  is 


^  =  PV'  (V  X  v* )  -  /*P[ V  X  (k  X  V* )]  +  P( V  X  7*k)  -  Re[V  x  (v*  •  Vv* )] 


-[V  X  (v*  .  Vv* )]  -  [V  X  (v*  •  Vv* )]  -  Re 


Vx  w' 


av* 

dz 


(2.33) 


We  specify  appropriate  boundary  conditions  for  the  system  of  equations  (2.24)  and 
(2.33).  For  simplicity,  the  lateral  boundaries  are  cyclic,  and  the  upper  vertical  boundary  is 
rigid,  stress-free,  and  perfectly  conducting.  In  contrast,  an  important  source  of  energy  for 
the  circulations  exists  at  the  lower  boundary;  these  boundary  conditions  are  specified 


using  surface  layer  similarity  theory  as  described  in  Zuccarello  (1994).  The  vertical 
boundary  conditions  are  summarized  below  for  5^.  >  0  and  s^>0,  which  is  consistent 

with  the  flow  depicted  in  Fig.  2.1: 
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dz 

(2.34) 

(2.35) 

(2.36) 

o 
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(2.37) 

z  dz 

(2.38) 

where  ,  the  Schramm  momentum  constant,  and  Sj ,  the  Schramm  temperature  constant, 
are  the  similarity  parameters  as  defined  in  Zuccarello  (1994).  These  parameters  control 
the  forcing  at  the  lower  vertical  boundary.  For  consistency,  the  mean  wind  profile  at 
z*  =0  should  satisfy  (2.26),  (2.35)  and  (2.36).  As  a  consequence  of  (2.26),  (2.35)  and 
(2.36),  we  get 


9w*(0)  3^w*(0) 

“a?  dz^ 


=  0 


(2.39) 


Note  that  this  equation  does  not  handle  the  vertical  velocity  at  the  lower  boundary  as  we 
would  like,  but  rather  it  involves  derivatives  of  the  vertical  velocity. 
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2.1.2.  The  Spectral  System 

Zuccarello  used  the  Galerkin  technique  or  spectral  method  to  form  a  spectral 
model.  We  use  Fourier  expansions  for  the  variables  that  were  determined  by  a  suitable 
eigenvalue  problem  using  the  above  boundary  conditions.  The  Galerkin  method  applied 
to  (2.24)  and  (2.33)  yields  a  nonlinear  dynamical  system  (NDS).  We  then  truncate  the 
NDS  to  60  equations  based  on  20  variable  amplitudes  for  temperature  and  40  for  velocity. 
Each  partial  differential  equation  term  may  be  represented  by  the  product  of  a  60x  60 
matrix  and  a  60-element  vector  representing  the  60  amplitude  coefficients.  The  resulting 
NDS  takes  the  form 


Ay  =  By  +  C(y)y,  (2.40) 

where  A  is  a  60x  60  invertible  matrix  of  inner  products  multiplying  the  60-element 
vector  y  that  contains  the  temporal  derivatives  of  the  amplitude  coefficients,  5  is  a 

60x  60  constant  matrix,  arising  from  the  linear  terms  of  the  partial  differential  equations, 
multiplying  the  60-element  vector  y  that  contains  the  amplitude  coefficients,  and  C  is  a 

60x  60  linear-in-y  matrix,  based  on  the  advective  terms  of  the  partial  differential 
equations,  multiplying  the  same  60-element  vector  y.  Because  A  is  invertible,  we  may 
solve  for  the  time-dependent  amplitude  coefficients  by  numerically  integrating  the 
following  system: 


y  =  A-'By+A-'Ciy)y. 


(2.41) 
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The  truncated  Fourier  expansion  for  the  dependent  variable  T*  is  composed  of 
temporally  dependent  amplitudes,  spatially  dependent  trigonometric  functions  with 
horizontal  wavenumbers  I  and  m ,  and  vertical  basis  functions: 

T\x\y*,z*,t*)  =  '^'^Ty  {t*)trig.  {x\y*)Fj(z*),  (2.42) 

1=0  j=i 

where 


(2.43) 

fn'g,  (x* ,  y * ]  =  sin  Inlx*  sin  lumy* 

(2.44) 

trig^ix*  ,y*^  =  sm2%lx*  coslnmy* 

(2.45) 

trig^  (x* ,  y * )  =  cos  2%lx*  sin  2Kmy* 

(2.46) 

trig^i^x*  ,y*^  =  cos27t/x*  cos27tmy‘ 

(2.47) 

F,{z*)  =  coshO),z’ - ^ — sinhcOjZ*  for  <  1 

5j-C0, 

(2.48) 

fj(z*)  =  cosco,z* - ^ — sin®,z*  for  >  1 

(2.49) 

57.(0,  =  tanho),  for  Sj  <  1 

(2.50) 

57-CD,  =  tanco,  for  Sj  >  1 

(2.51) 

F(z*)  =  cosa),z* - ^ — sin®,z*  for  j  =  2,3,4 

Sjdij 

(2.52) 

Sj(£)  j  =  tan  CO  ^  for  j  =  2,3,4, 

(2.53) 

and  the  vertical  basis  functions  are  orthogonal. 

The  expansion  for  the  velocity  vector  v*  is  more  complex.  Because  it  is 
nondivergent,  v*  is  represented  (nonuniquely)  by  the  curl  of  a  vector  streamfunction  A* 
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that  is  the  sum  of  two  vector  streamfunctions  V|/*  andT|* .  Therefore,  we  express  the 
velocity  perturbation  vector  v*  as 

V*  =  VxA*  =  Vx\|/*  +  Vxri*.  (2.54) 

We  use  a  truncated  Fourier  expansion,  similar  to  that  for  T* ,  to  represent  the  dependent 
vector  streamfunctions  \|/*  andrj* : 

\|/*(x*,/,z‘,t‘)  =  +  +  (2.55) 

p=0 

ri‘(x*,/,z‘,t*)  =0i  +  2^^rip^(t*)rng^(j:‘,/)/i^(z*)j  +  0k,  (2.56) 

p=0  q=] 

where  is  defined  in  (2.43)-(2.47),  and 

/i^(z*)  =  cosG3^z*-j„,C3^sin03^z*  for  ^  =  1,2, 3, 4  (2.57) 

SJJ5  com  ^  for  ^  =  1,2, 3, 4,  (2.58) 

where  the  vertical  basis  functions  h^[z*)  are  nonorthogonal. 

During  the  Galerkin  process,  the  expansions  are  put  into  the  appropriate  terms  in 
the  dimensionless  partial  differential  equations,  multiplied  by  appropriate  individual  basis 
functions,  and  then  simplified  by  integrating  the  trigonometric  functions  in  the 
expansions  and  using  appropriate  identities. 
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2.2.  Energetics  Analysis 


In  Zuccarello  (1994),  it  is  hypothesized  that  computer  runs  of  the  model  there  fail  to 
converge  because  of  an  unknown  energy  source.  Here,  we  performed  an  energetics 
analysis  to  determine  all  the  energy  sources  and  sinks  in  the  system.  Using  the 
thermodynamic  equation  and  the  equation  of  motion,  we  can  determine  the  energy  rate 
equations  for  the  system.  These  energy  rate  equations  are  composed  of  obvious  energy 
sources  and  sinks,  as  well  as  boundary  terms  that  can  be  either  sources  or  sinks.  By 
properly  choosing  the  boundary  conditions,  we  can  eliminate  most  of  the  boundary  terms 
and  so  simplify  the  energy  rate  equations.  Using  energetics  analysis  as  a  diagnostic  tool, 
we  determine  whether  our  current  boundary  conditions  adequately  control  the  system 
energetics.  Appendix  A  fully  describes  the  energetics  analysis. 

In  the  first  step  in  the  energetics  analysis,  we  calculate  the  kinetic  energy  {KE)  rate 
equation  from  the  equation  of  motion.  We  multiply  the  equation  of  motion  (2.25)  by  the 
velocity  vector  v*  and  then  take  the  volume  integral  of  the  product.  The  result  of  the 
vector  dot  product,  integration  by  parts,  introduction  of  the  boundary  conditions  (2.34)  - 
(2.39),  and  integration  over  z  when  possible,  is: 


B 


(1)  (2)  (3) 


PT*w*dV+ 


V 


(2.59) 


(4) 


(5) 


(6) 
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where  KE  is  the  KE  rate  of  change  and  B  is  the  horizontal  boundary. 

The  six  terms  on  the  right  side  of  (2.59)  represent  sources  and  sinks  of  energy.  The 
first  term  is  a  boundary  term  due  to  pressure  work.  The  second  term  is  traditionally  the 
KE  dissipation  term  and  is  clearly  an  energy  sink.  The  third  and  fifth  terms  are  also 
boundary  terms.  The  three  boundary  terms  (1),  (3),  and  (5)  are  either  energy  sources  or 
sinks  depending  on  the  sign  of  the  velocity  term  within  the  integral  and  the  sign  of  the 
Schramm  constant  in  term  (3)  (for  this  model  >  0).  The  fourth  term  is  a  convective 
energy  source  term  describing  the  conversion  of  AE  to  KE.  We  denote  it  as  C{AE,  KE) . 

The  last  term  (6)  is  a  KE  generation  or  Reynolds  stress  term  representing  the  dynamic 
source  of  energy. 

Next,  we  calculate  the  available  energy  (AE)  rate  equation  ( AE )  from  the 
thermodynamic  equation  (2.24)  and  the  C{AE,KE)  term  in  (2.59).  We  need  the  same 

conversion  term,  but  with  the  opposite  sign,  in  the  AE  equation.  Therefore,  we  add  and 

• 

subtract  the  required  term  in  order  to  introduce  it  appropriately.  We  derive  the  AE 
equation  by  multiplying  the  volume  integral  of  the  thermodynamic  equation  by  an 
appropriate  factor  to  obtain  —C(AE,KE)  from  the  modified  convection  term.  After 

integrating  by  parts,  applying  the  boundary  conditions  (2.34)  -  (2.39),  and  integrating  over 
z*  when  possible,  we  get  the  AE  equation: 


(2.60) 


(4) 


(5) 
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The  six  terms  on  the  right  side  of  (2.60)  represent  sources  and  sinks  of  AE.  The 
first  term  is  a  boundary  term  that  depends  on  the  sign  of  the  temperature  Schramm 
constant,  Sj .  We  use  5^  >  0  in  this  model,  so  term  (1)  is  an  AE  source.  The  second  term 

is  an  AE  dissipation  term  and  is  clearly  a  sink.  The  third  term  comes  from  the  original 
Rayleigh  term  in  the  thermodynamic  equation  and  our  correction  term.  Assuming  that 
r*  and  w*  should  be  positively  correlated,  the  third  term  is  an  AE  source  for  Ra>-1 

and  a  sink  for  Ra<-\.  The  case  of  Ra  =  -\  corresponds  to  a  thermodynamically 
neutral  atmosphere  as  no  thermal  energy  arising  from  the  mean  lapse  rate  A^T/zj.  can  be 
converted  to  KE.  In  term  (4),  we  find  the  C(AE,KE)  term  of  the  KE  equation,  but  with 

the  opposite  sign.  It  represents  a  source  or  sink  of  AE,  but  not  a  net  source  of  total  energy 
KE+AE.  The  last  term  is  another  boundary  term,  and,  depending  on  the  sign  of  the 
vertical  velocity  at  z*  =  0 ,  is  either  a  source  or  sink  of  AE. 

We  have  completed  the  energetics  analysis  for  the  original  equations  of  the 
problem,  but  the  model  uses  the  vorticity  equation  (2.33).  Therefore,  we  need  to  continue 
the  analysis  by  examining  the  KE  equation  derived  from  the  vorticity  equation.  For 
physical  consistency,  we  must  obtain  the  same  KE  equation  found  from  the  equation  of 
motion  (2.59).  However,  we  cannot  get  this  form  by  multiplying  (2.33)  by  the  velocity 
perturbation,  v* .  To  get  the  correct  form,  we  multiply  the  vorticity  equation  by  the 
vector  streamfunction  A*  and  then  take  the  volume  integral  (see  Appendix  A  for  details). 
The  result,  after  vector  multiplication,  integration  by  parts,  application  of  the  boundary 
conditions  (2.34)  -  (2.39),  and  integration  over  z*  when  possible,  is: 
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(2.61) 


where  the  tensor  A^^  is  the  streamfunction  component  V|/ *  when  k  =  \,  rj*  when  k  =  2, 
and  0  when  k  =  3. 

A  quick  comparison  of  (2.59)  and  (2.61)  shows  the  two  KE  equations  are  quite 
different.  However,  a  closer  look  reveals  some  interesting  results.  We  do  not  find  term 
(1),  the  pressure  boundary  term  of  (2.59),  directly  in  (2.61).  However,  we  find  the  next 
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five  terms  in  both  (2.59)  and  (2.61).  The  terms  are:  the  dissipation  term  (2);  a  boundary 
term  (3);  the  C[AE,KE)  term  (4);  another  boundary  term  (5);  and  the  KE  generation 

term  (6).  We  can  find  no  direct  comparison  in  (2.59)  with  the  six  remaining  terms  in 
(2.60).  Therefore,  these  six  terms  must  come  from  the  pressure  boundary  term  in  (2.58), 
since  these  two  equations  must  be  equal  and  the  pressure  boundary  term  is  the  only  term 
not  found  in  (2.61). 

The  last  six  terms  in  (2.61)  are  boundary  terms,  as  is  the  pressure  work  term  in 
(2.59).  As  stated  earlier,  if  we  choose  our  boundary  conditions  correctly  then  we  can 
eliminate  most  of  the  boundary  terms  in  the  energy  rate  equations.  Besides  these  related 
terms,  terms  (3)  and  (5)  are  boundary  terms  in  the  KE  equations,  and  term  (5)  is  a 
boundary  term  in  the  AE  equation.  Examining  all  the  boundary  terms  reveals  a 
dependence  on  the  vertical  velocity  w*  and  the  stream  function  Al  at  z*  =  0 .  We  chose 

the  simplest  new  boundary  condition  that  eliminates  most  of  these  terms: 


A;(z‘  =  0)  =  0.  (2.62) 

As  a  consequence  of  (2.62),  we  get 

w*(z‘=0)  =  0.  (2.63) 

Before  accepting  these  simple  and  obvious  boundary  conditions,  we  make  sure  that 
they  represent  observed  atmospheric  conditions.  It  could  be  argued  that  the  vertical 
velocity  perturbation  at  10  m  above  the  sea  surface  is  not  identically  zero.  However,  with 
typical  scaling  arguments,  we  can  show  that  the  vertical  velocity  perturbation  is  negligible 
when  compared  with  the  horizontal  velocity  perturbations.  Additionally,  we  do  not  need 
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a  nonzero  vertical  velocity  to  transfer  the  momentum  and  heat  flux  to  the  surface,  since 
smaller  scale  eddies  are  responsible  for  this  transfer. 

The  two  additional  boundary  conditions  help  solve  Zuccarello’s  problems  of  an 
unknown  energy  source  and  a  maximum  vertical  velocity  at  the  lower  boundary.  The 
new  boundary  conditions  eliminate  the  many  possible,  unrealistic  boundary  energy 
sources  and  they  control  the  vertical  velocity  at  the  lower  boundary.  The  first  boundary 
condition  (2.62)  also  helps  to  move  the  maximum  momentum  flux  higher  into  the 
boundary  layer.  The  horizontal  velocities  are  curls  of  the  two  streamfunctions  in  (2.54) 
and  have  z-derivatives  of  the  streamfunction  in  their  definitions.  Therefore,  boundary 
condition  (2.62)  does  not  force  the  horizontal  velocities  to  zero,  but  it  does  reduce  them. 
This  reduction  of  horizontal  velocity  and  our  control  of  the  vertical  velocity  at  the  lower 
boundary  forces  the  maximum  momentum  flux  to  occur  higher  in  the  boundary  layer. 

The  addition  of  these  two  boundary  conditions  reduces  the  energetics  to; 
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The  terms  in  the  KE  equation  are:  KE  dissipation  (1),  a  boundary  KE  sink  term  (2)  (since 
>  0),  the  C[AE,KE)  term  (3),  and  the  KE  generation  term  (4).  Similar  terms  are  in 

the  AE  equation:  a  boundary  AE  source  term  (1)  (since  >  0 ),  the  AE  dissipation  term 
(2),  the  AE  generation  term  (3),  and  the  C[AE,  KE)  term  (4).  Note  that  boundary  layer 

flow  driven  by  heat  flux  from  the  bottom  and  modified  only  by  shear  within  the  boundary 
layer  occurs  when  Ra  =  -l.  We  examine  this  special  case  in  Chapter  3. 

2.3.  Incorporation  of  New  Boundary  Conditions 


The  addition  of  two  new  boundary  conditions  (2.62)  -  (2.63)  to  the  original  Monin- 
Obukhov  boundary  conditions  (2.34)  -  (2.39),  greatly  changes  our  system  of  equations. 
However,  we  did  not  want  to  completely  rederive  and  rewrite  the  model  at  this  time. 
Instead  we  were  able  to  adjust  the  existing  system  with  some  minor  changes  to  the 
vertical  expansion  functions  for  the  streamfunctions  that  now  satisfy  both  new  boundary 
conditions. 

We  apply  the  new  boundary  conditions  at  the  lower  boundary  z*  =  0 .  So  we  start 
by  examining  the  vector  streamfunction  A*  using  (2.55)  and  (2.56): 


A*  (x‘ ,  y  ‘ ,  ,  r*)  =  ^  [{\|r  (x‘ ,  y 

p=0  q=\ 


)}j+0k]. 


(2.66) 
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The  only  dependence  on  z*  in  (2.66)  is  in  the  vertical  basis  function  h^[z*) .  If  we  want 
A*(z*  =  O)  =  0 ,  then  we  must  force  h^{0)  =  0 ,  for  each  of  the  four  vertical  wavenumbers 

q- 

However,  using  the  definition  of  h^{z)  in  (2.57),  we  find  its  value  for  each  of  the 
four  vertical  wavenumbers  at  the  lower  boundary  is: 


(2.67) 


From  (2.67)  we  see  a  very  simple  way  to  make  the  new  vertical  basis  functions  vanish 
leading  to  A*(z*  =  O)  =  0 .  We  define  new  vertical  basis  functions  ,  for  q'  =  1,2,3 , 

by  taking  appropriate  differences  of  the  old  basis  functions  ; 


h 

new] 

h 

^%ew2 

Kewi 


“  ^olc/l  ”■  Ki(12 

(2.68) 

~  Kidi  ~  Kid3 

(2.69) 

~  Kld?>  Kid  A  ‘ 

(2.70) 

The  three  new  vertical  basis  functions  reduce  the  number  of  vertical  wavenumbers  in  the 

streamfunction  representation  from  four  to  three.  We  use  (2.57)  to  define  the  old  basis 

functions  and  then  use  these  values  to  obtain  the  new  basis  functions  (2.68)  -  (2.70). 

Once  put  into  a  new  Galerkin  expansion,  the  resulting  new  function  satisfies  the  boundary 
condition  A*(z‘  =  O)  =  0 : 


+{'n  pci’  ‘  y^'^gp  Kc.,-  (2* )} j + Ok  . 


(2.71) 
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Substituting  (2.68)-(2.70)  in  (2.71)  yields 


p=Oq=] 


P‘1'  (^*  (x* ,  y){Kw  (^* )  -  (^* ))} j + Ok  . 


(2.72) 


At  z*  =  0 ,  the  old  basis  functions  =  1  by  (2.67)  and  therefore  their  differences  in 
(2.72)  cancel.  Summing  over  p  and  q'  results  in  A*(z*  =  O)  =  0 . 

When  we  calculate  the  curl  of  A*  in  (2.71),  we  get  the  new  velocity  component 
equations.  The  vertical  velocity  expansion  is: 


w 


^=0^'=lL  OX 


(2.73) 


Summing  over  p  and  q' ,  we  get  w*(0)  =  0 ,  which  is  our  other  new  boundary  condition. 

Appendix  B  offers  a  closer  look  at  incorporating  these  new  boundary  conditions  into  our 
system  of  equations. 


2.4.  The  Computer  Model 


The  computer  code  used  to  represent  the  physical  and  mathematical  model 
presented  in  the  previous  sections  is  written  in  the  FORTRAN  programming  language.  It 
is  compiled  using  the  WATFOR77  compiler.  The  first  part  of  the  code  determines  the 
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values  for  co  j  for  j=  1,2,  3, 4,  and  03  ^  for  ^  =  1 ,  2,  3, 4,  based  on  input  values  of  the 
constants  5^  and  Sj. .  The  values  of  these  constants  are  provided  by  the  subroutine 
discussed  in  Appendix  A  of  Zuccarello  (1994).  An  iterative  technique  is  used  to  solve 
(2.50)-(2.51),  (2.53),  and  (2.58)  for  the  eigenvalues.  Zuccarello  suggests  that  for  new 
values  of  and  s-j. ,  a  new  set  of  eigenvalues  be  put  into  the  model  as  a  seed  set  for  the 

iterative  process.  However,  after  running  several  different  cases  with  different  values  of 
and  Sj ,  we  find  the  iterative  process  adequately  yields  the  new  eigenvalues  from  a 

single  seed  set.  Prior  to  integrating  the  system,  we  choose  values  for  the  parameters 
l,m,a^,a^.,Ra,  Re,  andP. 

The  above  parameter  values  are  used  to  fill  the  matrices  representing  the  terms  in 
the  NDS.  The  individual  matrices  that  compose  A,  B,  and  C  in  (2.39)  are  grouped  as 
temporal  derivative  term  matrices,  linear  term  matrices,  and  nonlinear  term  matrices, 
respectively.  The  temporal  derivative  term  matrix  A  and  the  linear  term  matrix  B  are 
determined  prior  to  the  integration  because  A  and  B  do  not  depend  on  the  spectral 
coefficient  in  y .  An  IMSL  routine  computes  A”'  and  a  different  IMSL  routine  computes 

.  The  nonlinear  term  matrix  must  be  calculated  during  each  integration  time  step 
because  C  depends  on  y .  The  calculation  of  A“'C(y)  may  also  be  done  during  each  time 

step. 

During  a  preliminary  linear  analysis  of  the  model,  we  discovered  a  problem  with 

.  There  were  some  terms  in  the  matrix  that  should  be  zero  because  of  algebraic 
cancellation.  However,  owing  to  round-off  errors,  the  values  of  these  terms  had  values  on 
the  order  of  10”'®  to  10"'® .  In  scaling  arguments,  we  would  normally  take  these  numbers 
to  be  zero,  but  the  IMSL  subroutine  we  use  for  the  linear  analysis  did  not.  When  we 
perform  eigenvalue  analyses  of  the  temperature  and  vorticity  dissipation  terms,  which 
appear  only  on  the  diagonal  of  the  matrix,  we  expected  to  get  only  real  eigenvalues. 
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Instead  we  obtained  complex  eigenvalues  because  of  the  small  nonzero  terms  in  the 
matrix.  We  solved  this  problem  by  copying  all  the  terms  that  were  supposed  to  be 

nonzero  into  another  matrix  filled  with  zeros. 

When  we  defined  the  new  vertical  basis  function  ,  we  reduced  the 

streamfunction  vertical  wavenumbers  from  four  to  three.  This  result  reduces  the  number 
of  model  equations  and  coefficients  from  60  to  50.  Therefore,  we  rewrote  the  parts  of  the 
model,  involving  the  vertical  basis  functions,  using  the  new  definitions  (2.68)  -  (2.70)  and 
construct  A ,  B  and  C  as  50  x  50  matrices. 

The  IMSL  integration  routine  uses  the  Adams-Gear  method  and  requires  user- 
supplied  initial  conditions  for  the  50-element  vector  y .  At  each  time  step,  an  IMSL 
routine  calculates  A~'By  and  A'’C(y)>'.  The  resulting  50-element  vectors  from  the 

linear  and  nonlinear  matrices  are  added  together  to  give  the  right  sides  of  the  50  equations 
for  which  we  are  solving. 

The  integration  routine  uses  these  results  to  calculate  new  values  for  the  50  time- 
dependent  amplitude  coefficients  after  each  time  step.  We  integrate  the  model  until  the 
square  root  of  the  sum  of  the  squares  of  the  coefficients  equilibrates.  After  the  model 
integration  is  complete,  these  amplitude  coefficients  are  used  to  construct  the  model 
solution  using  (2.42),  (2.54)  and  (2.71).  Examples  of  this  solution  are  given  in  the 
following  chapter. 
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Chapter  3 

HI-RES  CASE  STUDY 


We  evaluate  the  model  with  the  new  boundary  conditions  by  performing  a  series 
of  model  analyses  and  integrations  using  observed  marine  atmospheric  boundary  layer 
(MABL)  data  from  1993  Hi-Res  2  field  experiment  to  calculate  the  model  parameter 
values.  This  experiment  yielded  coincident  ship  measurements  of  the  atmospheric  and 
oceanic  state  and  synthetic  aperture  radar  (SAR)  imagery  from  ERS-1  (Sikora  et  al. 

1995).  The  analyses  include  linear  studies  using  different  physical  parameter  values  and 
wind  profiles.  We  determine  the  preferred  physical  parameter  values  for  our  case  study 
using  linear  eigensystem  evaluations  (Laufersweiler  1987).  The  examination  of  the 
preferred  physical  parameter  values  yields  a  set  of  aspect  ratios  associated  with  a 
bifurcating,  temporally  periodic  solution.  We  then  perform  some  integrations  using 
several  sets  of  these  preferred,  physical  parameter  values  and  initial  conditions.  After  the 
integrations,  the  analyses  include  energy  conservation  checks  on  the  model  solutions  and 
graphical  representations  of  the  evolution  of  the  solution.  We  use  these  integrations  to 
locate  code  errors  and  to  fine-tune  the  model  parameters. 

Using  one  set  of  aspect  ratios,  we  perform  the  Hi-Res  2  case  study  to  investigate 
the  solutions  for  boundary  layer  spanning  eddies  (ELSE).  These  eddies  are  composed  of 
component  circulations  parallel  to  the  mean  wind  and  perpendicular  to  the  mean  wind; 
both  roll  circulations  and  cellular  ones  are  represented.  We  examine  the  ELSE  surface- 
stress  patterns  produced  by  the  model  and  compare  them  with  the  pattern  observed  in  our 
case  study  Hi-Res  2  SAR  imagery.  We  also  compare  the  results  with  those  of  other 
models  to  assess  both  the  performance  of  our  model  and  the  effects  of  incorporating  the 
forcing  rates  in  the  lower  vertical  boundary  conditions.  We  describe  the  Hi-Res  2  data  in 
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Section  3.1  and  present  the  case  study  data  and  analyses  yielding  the  preferred  physical 
parameter  values  in  Section  3.2.  The  results  of  this  case  study  are  given  in  Section  3.3. 

3.1.  Description  of  Hi-Res  2  Data 

The  observed  data  for  our  case  study  are  from  the  Hi-Res  2  field  experiment 
conducted  during  June  1993.  The  Hi-Res  cruise  coincided  with  SAR  imaging  from  the 
European  Remote  Sensing  Satellite  ERS-1  off  the  coast  of  North  Carolina.  From  16  June 
to  18  June  1993,  the  RA^  Columbus  Iselin  and  the  RA^  Bartlett  collected  observations  at 
the  sea-surface  and  at  z  =  lOm  within  the  area  enclosed  by  the  coordinates  35.2°  -  36.6° 

N  and  73.8°  -  75.0°  W.  We  examined  an  ERS-1  SAR  composite  image,  published  by 
Sikora  et  al.  (Figure  1,  1995),  dated  1538  GMT  17  June  1993  that  covered  this  same  area. 

Sikora  et  al.  (1995)  found  two  distinct  patterns  on  this  image.  They  concluded 
that  the  mottled  pattern  on  this  image,  south  of  the  Gulf  Stream  North  Wall  (GSNW),  is 
the  signature  of  convectively  driven  ELSE,  while  the  marbled  pattern,  north  of  the 
GSNW,  is  the  signature  of  a  stable  MABL.  Using  these  criteria,  we  chose  to  study  the 
convective  region  near  35.8°  N,  74.2°  W  that  displays  an  alternating  dark/light  pattern 
with  an  apparent  kilometer-scale  wavelength.  The  R/Y  Columbus  Iselin  researchers 
recorded  surface  observations,  kindly  shared  with  us  by  Dr.  Jim  Edson  of  Woods  Hole 
Oceanographic  Institution,  in  this  region  near  the  time  of  the  SAR  image. 

The  meteorological  conditions  in  the  Hi-Res  region  at  this  time  were  relatively 
quiet.  According  to  Sikora  et  al.  (1995),  the  1500  GMT  17  June  1993  Nationzd 
Meteorological  Center  (NMC)  surface  analysis  indicates  that  the  imaged  area  was  located 
within  a  uniform,  relatively  warm  air  mass  south  of  a  stationary  front.  The  NMC  upper 
air  analysis  from  1200  GMT  17  June  1993  shows  the  imaged  area  to  be  just  east  of  a 
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large,  upper- air  ridge  (Sikora  et  al.  1995).  Also,  Sikora  et  al.  (1995)  report  that 
observation  notes  from  meteorological  personnel  aboard  the  RA^  Columbus  Iselin 
indicate  a  dramatic  atmospheric  meso-front  occurred  in  the  vicinity  of  the  GSNW. 

The  absence  of  any  synoptic-scale  disturbances  suggests  that  locally  driven  ELSE 
are  responsible  for  the  imaged  stress  patterns  (Sikora  et  al.  1995).  The  RA^  Columbus 
Iselin  data  from  over  the  Gulf  Stream  displays  a  distinct  air/sea  temperature  difference, 
with  sea-surface  temperatures  (SST)  greater  than  the  air  temperatures  in  the  region  of 
marbled  sea-surface  patterns.  This  air/sea  temperature  difference  is  one  of  the 
mechanisms  for  ELSE  described  in  Chapter  2  and  modeled  here. 

3.2.  Case  Study  Data  and  Preferred  Parameter  Values 

The  meteorological  staff  onboard  the  RA^  Columbus  Iselin  took  observations  an 
average  of  three  to  four  times  an  hour  on  17  June.  We  average  these  observations  for 
each  hour  (Table  3.1).  The  RfW  Columbus  Iselin  was  in  the  vicinity  of  our  case  study 
point,  35.8°  N,  74.2°  W,  from  1400  GMT  17  June  through  0000  GMT  18  June.  Since 
observations  were  taken  at  10  m,  we  chose  z  =  =  10  m .  The  observed  wind 

measurements  were  relative  to  the  fixed  Earth.  We  subtract  the  observed  current  from  the 
observed  winds  to  make  the  wind  velocity  relative  to  the  ocean  —  U  Rel  and  V  Rel.  We 
infer  the  mixing  ratio  at  the  sea  surface  by  assuming  that  the  air  is  saturated  at  this  level. 
We  use  the  averages  of  all  the  observed  parameter  values  over  each  hour  to  specify  the 
model  parameter  values  for  a  case.  Table  3.1.  (Note:  the  observed  data  shared  with  us  did 
not  include  any  1700  GMT  and  2000  GMT  17  June  data.) 


Table  3.  L  Hi-Res  2  Observed  MABL  data.  Note:  1700  GMT  and  2000  GMT  data  not  available. 
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The  physical  parameter  values  from  Table  3.1  —  U Rel,  V Rel,  SST,  T\0,  qsea, 
qlO  —  are  input  into  the  subroutine  described  in  Appendix  A  of  Zuccarello  (1994).  This 
subroutine  determines  values  for  the  surface  roughness,  Zq  ,  and  the  Monin-Obukhov 

length,  L.  The  values  obtained  for  Zq  and  L  are  then  used  in  another  subroutine  that 
calculates  values  for  the  two  constant  lower  boundary  Schramm  forcing  parameters, 
and  Sj .  The  values  for  Zq  ,  L,  and  Sj-  calculated  from  the  different  Hi-Res  data  points 

are  found  in  Table  3. 1 .  Our  values  of  L  agree  with  those  in  Sikora  et  al.  (1995).  The 
range  of  values  of  s„  and  Sj-  for  17  June  are  exponentially  related,  as  shown  in  Fig.  3.1. 

From  Table  3.1,  we  chose  the  averaged  1900  GMT  data  for  our  case  study. 

During  this  time  period,  the  RA^  Columbus  Iselin  was  in  the  middle  of  the  dark/light 
pattern  on  the  SAR  image.  This  pattern  is  oriented  35°  west  of  north  and  has  horizontal 
spacing  of  approximately  1000  m.  The  parameter  values  for  and  Sj-  (darkened  square 

in  Fig.  3.1)  are  typical  of  those  values  characterizing  the  state  of  the  17  June  MABL.  The 
relative  winds  at  z  =  10  m  were  2  ms'^  from  the  southwest.  We  plot  the  value  of  the 

relative  winds  and  the  air/sea  temperature  difference  on  Woodcock's  (1975)  regime 
diagram  in  Fig.  3.2.  The  regime  diagram  indicates  that  the  1400-2400  GMT  averaged 
data  are  in  the  three-dimensional  cell  region  (1900  GMT  point  is  darkened  diamond). 

The  1900  GMT  upper  air  sounding,  taken  from  the  RA^  Columbus  Iselin  during  this  time, 
showed  an  inversion  base  height  of  900  m.  This  inversion  base  height  defines  our  upper 
boundary  height  Zy-  =  900  m .  The  expected  aspect  ratios  for  the  observed  spacing  on  the 
SAR  image  are  thus  =  0.73  and  a,,  =  0.5 1 .  The  sounding  winds  were  not  available 

through  the  boundary  layer.  However,  researchers  on  board  the  IW  Columbus  Iselin 
noted  that  the  wind  speeds  were  light  through  the  boundary  layer. 


Current-Relative  Wind  Speed 
@  1 0  m  (m/s) 
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Rather  than  use  a  constant  wind  profile,  we  use  formulas  which  we  may  fit  to 
natural  boundary  conditions  and  to  winds  observed  at  the  surface  and  aloft.  Since  the 
winds  aloft  are  not  known,  we  use  the  following  formulas  for  the  profiles  Uiz)  and  V(z) 


U{z)  =  a,  z+by  z^'^+c^  (3.1) 

y(z)  =  ayZ+by  z"'^+Cy  z'^'  +dy .  (3.2) 


We  use  (2.14)  -  (2.23)  to  make  these  profiles  dimensionless.  We  require  that  the  mean 
wind  profiles  satisfy  the  same  boundary  conditions  that  the  perturbation  wind  velocity 
satisfies.  These  boundary  conditions  are 


oz  az 


U\z  =0)  =  URell  y*(z*  =  0)  =  VRel  /  |V, 

=  I)  =  Ur.,  /|V„  I,  V  (z-  =  1)  =  Vr„  /  |V, 


Top[ 


'Top\ 


(3.3) 

(3.4) 


Since  the  mean  wind  varies  only  in  the  vertical  direction,  it  is  also  incompressible.  We 

obtain  the  Reynolds  number  Re  through  its  definition  (2.27).  In  (3.3)  -  (3.4)  we  then  use 
the  values  from  Table  3.1  for  ,  U Rel,  and  V Rel  at  z  =  hi^  =  \Qm  (z*  =  O)  to  solve  for 

the  variables  a^j  ,bjj ,  c^, ,  ay  ,by ,  and  Cy  in  (3. 1)  and  (3.2).  The  variables  and  dy  are 
free  variables  that  we  use  to  control  the  shear  and  turning  in  the  profile.  We  use  an 
EXCEL  spreadsheet  to  calculate  and  plot  the  profiles  for  different  values  of  Re  and  the 
free  variables  d^  and  dy  prior  to  inputting  them  into  the  model. 


Table  3.2.  Parameter  values  used  in  case  study. 


PARAMETER 

VALUE 

HEIGHT  OF  TOP  OF  DOMAIN  (Zt) 

900  m 

HEIGHT  OF  BOTTOM  OF  DOMAIN  (Hib) 

10  m 

ASPECT  RATIO  IN  THE  X-DIRECTION  (oj 

0.56 

ASPECT  RATIO  IN  THE  Y-DIRECTION  (a, ) 

0.38 

REYNOLDS  NUMBER  (Re) 

250.0 

U  MEAN  WIND  SPEED  hi  Z  =  Z^ 

1 .46  ms'^ 

U  MEAN  WIND  SPEED  AT  Z  = 

1.99  mY* 

V  MEAN  WIND  SPEED  AT  Z  = 

1.30  ms^ 

V  MEAN  WIND  SPEED  AT  Z  = 

1.99  ms'^ 

U  MEAN  WIND  PROFILE  PARAMETER  (du) 

1.40 

V  MEAN  WIND  PROFILE  PARAMETER  (dv) 

1.20 

GRAVITY  (;?) 

9.8  ms'^ 

CORIOLIS  PARAMETER  (/) 

0.0  y' 

EDDY  VISCOSITY  (v) 

1  n^s'^ 

EDDY  CONDUCTIVITY  (k) 

10  mV^ 

PRANDTL  NUMBER  (P) 

0.1 

SEA  SURFACE  TEMPERATURE  (Too) 

300.5  K 

TEMPERATURE  AT  Z  = 

25.4  °C 

MIXING  RATIO  AT  SEA  SURFACE 

22.33 

MIXING  RATIO  AT  Z  = 

12.10 

RAYLEIGH  NUMBER  (Rd) 

-1.0 

SCHRAMM  VELOCITY  CONSTANT  (Sn,) 

0.38 

SCHRAMM  THERMAL  CONSTANT  (st) 

0.91 
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Because  we  expect  the  ELSE  solutions  to  have  time  scales  on  the  order  of  an  hour 
or  less,  we  neglect  the  Coriolis  term  in  the  model.  In  our  model,  we  may  specify  only  one 
wavenumber  in  each  of  the  x-  and  y-directions,  and  so  we  set  /  =  1  and  m  =  1 .  The 
values  chosen  for  K  and  the  Prandtl  number  P  determine  the  value  of  the  v  using 
(2.30).  In  this  case,  k  =  and  P  =  0.1 ,  which  gives  v  =  1  m^s~' . 

For  the  next  step  in  finding  the  physical  parameter  values,  we  determine  the 
preferred  values  of  the  aspect  ratios.  The  aspect  ratios  describe  the  horizontal  domain 
shape  by  (2.30)  and  (2.31).  In  past  studies  of  nonlinear  dynamical  system  (NDS),  the 
preferred  aspect  ratios  have  been  found  by  performing  eigensystem  analyses  for  different 
values  of  the  forcing  represented  by  Ra  or  Re  (e.g.,  Shirer  1986;  Laufersweiler  and  Shirer 
1989,  1995;  Haack  and  Shirer  1992).  The  preferred  aspect  ratio  values  in  those  studies 
correspond  to  the  pairs  of  aspect  ratio  values  here,  and  their  minimum  Ra  or  Re  values 
that  are  given  by  a  conjugate  pair  of  eigenvalues  having  zero  real  part  in  the  linear 
stability  analysis,  corresponds  to  the  minimum  Re  value  here  (Laufersweiler  1987).  This 
result  signals  a  Hopf  bifurcation  producing  a  periodic  solution  that,  if  supercritical  and 
hence  stable,  is  realized  after  the  system  transients  have  disappeared  {e.g.  Pyle  1987). 

Figure  3.3  is  a  typical  example  of  a  plot  used  to  determine  the  preferred  values  of 
the  aspect  ratios.  The  curve  represents  the  vanishing  of  the  greatest  real  part  among  the 
eigenvalues.  The  lowest  point  on  this  curve,  giving  the  minimum  critical  Ra  or  Re  and 
the  preferred  aspect  ratio,  can  be  found  approximately  by  varying  the  aspect  ratios  in  an 
eigensystem  analysis  for  fixed  values  of  Ra  and  Re.  The  aspect  ratios  that  produce  the 
largest  positive  real  part  in  the  eigenvalues  for  this  fixed  Ra  and  Re  correspond  to  a  point 
approximately  above  the  preferred  aspect  ratio.  As  Ra  or  Re  is  varied  to  make  the 
greatest  real  part  tend  towards  zero,  the  corresponding  aspect  ratio  approaches  the 
preferred  value. 
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Fig.  3.3.  Schematic  plot  of  Re  versus  one  aspect  ratio  while  the  other  aspect  ratio  is  fixed  and  /?a  =  -1 . 
The  lower  curve  represents  the  vanishing  of  the  real  part  of  the  eigenvalue  having  the  largest  real 
part.  The  real  part  of  this  last  eigenvalue  is  negative  below  the  curve  and  positive  above  it.  The 
straight  line  is  given  by  the  value  of  Re  used  in  a  linear  analysis.  The  approximate  preferred  aspect 
ratio  can  be  found  numerically  by  determining  the  largest  real  part  along  this  line.  This  point 
represents  the  greatest  distance  from  the  line  to  the  zero  curve  and  therefore  is  above  the  well.  As 
Re  is  varied  to  make  the  largest  real  part  tend  toward  zero,  the  approximate  preferred  aspect  ratio 
approaches  the  true  preferred  aspect  ratio.  A  second  pair  of  eigenvalues  has  positive  real  part 
above  the  curve  for  the  second  bifurcating  mode.  The  Re  =  constant  line  is  chosen  below  this 
second  curve. 
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However  in  our  NDS  we  have  separated  the  thermal  forcing  at  the  lower  boundary 
from  the  thermal  forcing  within  the  boundary  layer.  The  Ra  value  in  our  NDS  controls 
the  thermal  forcing  within  the  boundary  layer  and  Sj  controls  the  thermal  forcing  at  the 
lower  boundary.  The  observed  meteorological  conditions  determine  the  value  of  Sj .  We 
choose  a  thermally  neutral  atmosphere  for  this  case  study,  so  Ra  =  -1  according  to  the 
energetics  analysis  discussed  in  Chapter  2  and  Appendix  A. 

Since  we  have  constant  Sj,  s^,  and  Ra,  but  an  unknown  wind  profile,  we 

determine  the  preferred  values  of  the  aspect  ratios  by  using  a  critical  Re,  instead  of  a 

critical  Ra,  in  our  eigensystem  analyses.  For  a  given  Re,  we  perform  a  linear  eigensystem 
zinalysis  as  we  loop  through  the  aspect  ratios  —  a^  and  a^. .  Since  the  Hi-Res  case  study 

is  a  light-wind  case,  we  choose  an  appropriate  initial  Reynolds  Number  Re  =  250 ;  this 
corresponds  to  wind  speeds,  relative  to  the  ocean  current,  of  1.95  mh  at  z^h  =  10  m  to 
2.8 1  mis  at  Zj  =  900  m .  Figure  3.4  is  a  plot  of  the  results  for  this  linear  analysis.  The 
plot  shows  the  real  part  of  the  last  eigenvalue  (eigenvalues  ordered  with  respect  to  the  real 
part)  from  the  analysis  for  each  aspect  ratio  set,  which  as  we  note  above  gives  the 
preferred  values.  From  the  plots,  we  look  for  the  largest  positive  values  contoured  toward 
a  ringed  maximum.  Figure  3.4  displays  an  obvious  contoured  maximum.  We  perform 
another  linear  analysis  around  this  maximum  using  a  finer  resolution  loop  of  the  aspect 
ratios.  From  this  refined  analysis,  we  determine  the  aspect  ratios  at  the  maximum  point. 
For  this  case  study,  we  then  look  at  a  similar  plot  for  the  next  pair  of  eigenvalues  up  in 
the  list  to  ensure  that  there  is  not  another  bifurcating  mode  indicated  for  the  above  set  of 
preferred  aspect  ratios  and  Re  value.  If  there  is  another  bifurcating  solution,  then  we 
decrease  the  Re  until  we  isolate  a  set  of  preferred  aspect  ratios  with  only  one  bifurcating 
solution  (Fig.  3.3).  The  final  set  of  aspect  ratios  for  the  smallest  maximum  real  parts 
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together  with  the  Re  value  are  chosen  for  use  in  our  case  study.  The  analyses  yielded  a 
Reynolds  number:  Re  =  250 ,  preferred  aspect  ratios:  =  056  and  =  0.38  and  only 

one  bifurcating  solution.  Now  that  we  found  a  Re  value,  we  can  plot  the  mean  current- 
relative  wind  profile,  Fig.  3.5. 

We  determine  the  initial  conditions  by  using  the  eigensystem  analysis.  We  find 
the  real  (or  imaginary)  parts  of  the  eigenvector  associated  with  the  preferred  aspect  ratios 
and  eigenvalue  found  above.  These  values  tell  us  the  approximate  relative  magnitudes  of 
the  coefficients  in  the  nonlinear  branching  solution. 

As  noted  in  Chapter  2,  the  IMSL  integration  routine  uses  the  Adams-Gear 
method.  We  perform  a  series  of  model  integrations  to  determine  the  proper  value  for  the 
time  step.  The  time  step  is  reduced  until  the  number  of  internal  subroutine  calls  is  small 
enough  to  allow  accurate  solutions  at  each  step.  We  find  the  largest  possible  time  step  to 
be  1  second.  Next,  we  run  the  model  integration  for  a  number  of  hours,  to  allow 
sufficient  time  for  the  system  to  equilibrate.  We  determine  when  the  system  reaches  a 
quasi-equilibrated  state  by  plotting  the  square  root  of  the  sum  of  the  squares  of  the  50 
time-dependent  amplitude  coefficients. 

Starting  with  the  original  eigenvector,  the  model  integration  took  1 8  hours 
(32,400  time  units)  to  equilibrate.  A  plot  of  the  square  root  of  the  sum  of  the  squares  is 
given  in  Fig.  3.6.  If  we  had  multiplied  the  eigenvector  by  a  constant  to  obtain  values  on 
the  scale  of  the  model  amplitude  coefficients,  and  then  had  used  the  resulting  values  as 
initial  conditions,  then  we  believe  that  the  model  might  have  equilibrated  faster.  The 
constancy  of  the  energy  after  equilibration  matches  results  acquired  by  Shirer  (1986),  Pyle 
(1987),  and  Haack  and  Shirer  (1992). 

We  also  check  the  validity  of  the  model  by  examining  the  periodicity  of  the 
solution.  Using  a  linear  eigensystem  analysis,  we  determine  the  values  for  the  imaginary 
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parts  Im  of  the  eigenvalues.  We  calculate  the  limiting  period  of  the  bifurcating  solution, 
using  the  Im  for  the  last  eigenvalue,  by  (Pyle  1987), 

Period  =  (3.5) 

Im 


This  value  must  then  be  converted  to  dimensional  units  using  (2.19).  Based  on  our 
analysis,  Im  =  229.5.  Using  the  above  equation,  we  find  a  dimensionless  period  of 
2.738  X 10”^ ,  or  36.22  minutes. 

Given  our  small-speed-shear  wind  profile,  this  period  should  correspond  to  the 

speed  of  movement  of  the  cells  caused  by  advection  by  the  mean  wind.  To  find  this 
speed,  we  first  note  that  the  ratio  L^jL^  =  a^ja^.  serves  as  a  measure  of  the  orientation 

angle  6  of  a  ELSE  axis,  for  which  tan6  =  ;  here  6  is  the  angle  between  the  ELSE 

axis  and  east.  At  any  level,  the  background  wind  angle  <])  is  given  by 
tan(|)  =  y/C/  -V*lu* .  Thus  it  follows  that  the  angle  ([)  -0  between  the  background 
wind  vector  and  the  ELSE  axis  obeys  (Zuccarello  1994) 


tan((l)-6)  = 


tan  (j)  -  tan  0 
l-l-tan(j)tan0 


ay -a^U* 

ay ay 


(3.6) 


For  the  case  discussed  here,  a^  =  056  and  a^.  =  0.38 ,  angle  0  of  the  ELSE  axis  is  56° 

(north  of  east),  and  so,  with  the  layer  mean  values  U  =  1.98  mi"' ,  V  =  1.70  mi"' ,  and 
[v|  =  2.61  mi"' ,  the  angle  (j)  -0  between  the  wind  direction  and  the  BLSE-axis  is 

approximately  14° .  The  component  U  along  the  ELSE  axis  is  given  by 


tJ  =  U  COS0  +V  sin0  , 


(3.7) 
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and  y  across  this  axis  is  given  by 

y  =-i/sine+ycos0  .  (3.8) 

These  components  introduce  advection  periods  equal  to  L/V  and  L/U  ,  where  L  is  the 
ELSE  axis  L  =  {zj-  hj^)j -yjal  +  .  Here  L  is  approximately  1330  m;  which  produces 

periods  of  35.1 1  minutes  and  8.8  minutes,  respectively.  The  first  period  is  close  to  the 
value  of  36.22  minutes  given  by  (3.5).  These  periods  are  under  an  hour  and  so  are 
consistent  with  our  neglecting  of  the  Coriolis  term. 

We  next  check  these  estimates  for  the  periods  against  the  nonlinear  results  for  one 
of  the  time-dependent  amplitude  coefficients  7^,  given  by  the  numerical  integration.  The 

periods  seen  in  the  graph  in  Figs.  3.7  are  approximately  36  minutes  and  9  minutes, 
therefore  demonstrating  that  the  model  integration  is  accurately  reproducing  the  expected 
periodicities  of  the  solution.  The  other  ELSE  coefficients  also  have  periods  of 
approximately  36  minutes  and  9  minutes,  while  the  ELSE  modification  terms  become 
steady.  These  results  are  also  consistent  with  the  results  of  Shirer  (1986),  Pyle  (1987), 
and  Haack  and  Shirer  (1992). 

3.3.  ELSE  Patterns 

We  use  the  physical  parameter  values  in  Section  3.2  to  examine  the  chosen  case  of 
ELSE  for  the  time  of  the  SAR  image.  Our  goal  is  to  determine  whether  the  model  might 
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Fig.  3.7.  Time  series  of  one  of  the  time-dependent  temperature  amplitude  coefficients  (TJ, ), 

demonstrating  that  the  bifurcating  solution  has  periods  of  approximately  36  minutes  and  9  minutes. 
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explain  the  patterns  seen  in  the  SAR  image  and  to  verify  that  the  new  boundary 
conditions  produce  the  expected  behavior  in  the  various  physical  quantities.  We  analyze 
and  compare  profiles  and  cross-sections  with  results  from  other  models,  and  then  we 
compare  planview  patterns  from  the  model  with  those  we  find  on  the  SAR  image. 

Our  first  step  is  to  determine  whether  the  model  captures  the  proper  spatial 
organization.  One  way  of  examining  this  structure  is  to  analyze  the  vertical  cross-section 
of  the  domain  along  the  two  axes:  x*  -along,  given  by 


X*  =  X*  COS0  -1- y*  sin 9  , 

(3.9) 

and  y'*  -across,  given  by 

y*  =  -X*  sin0  -l-  y*  cos0  . 

(3.10) 

These  axes  depend  on  the  orientation  angle  0  =  tan”'(a^/aj,) ,  and 

they  are  approximately 

along  and  perpendicular  to  the  mean  wind  direction,  respectively. 

Then  we  rotate  the 

wind  components  and  streamfunctions  via 

u*  =  m*cos9  +v*  sin9  . 

(3.11) 

V  *  =  —u  sin0  +  V*  COS0  . 

(3.12) 

xj/*  =\|/*cos9-r|*sin0 

(3.13) 

fi*  =\(/*sin0  +ti*cos9  . 

(3.14) 

This  rotation  does  produce  cross-sections  exhibiting  simple  horizontal  periodicity.  We 
present  the  contoured  streamfunctions  xjt*  and  fi*  in  the  cross-wind  direction  y*  in  Fig. 
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3.8(a)  and  Fig.  3.8(b),  and  in  the  along-wind  direction  x*  in  Fig.  3.9(a)  and  Fig.  3.9(b). 
We  immediately  see  that  the  ELSE  are  three-dimensional  as  predicted  from  the  data  on 
the  Woodcock  diagram  (Fig.  3.2).  In  addition,  the  circulations  fill  the  entire  vertical 
domain.  It  is  thus  clear  that  vertical  wavenumber  one  is  dominant,  but  there  is  also 
evidence  of  the  other  vertical  wavenumbers.  The  horizontal  trigonometric  functions  are 
interacting,  as  seen  by  the  lateral  nonsymmetry.  As  noted  above,  the  approximately 
repeating  horizontal  pattern  is  due  to  our  axis  rotation  by  the  orientation  angle  0  .  In  the 
tjl*  diagrams  (part  (a)  in  both  figures),  we  observe  large  horizontal  components  along  the 

lower  boundary  that  contribute  to  the  downward  transfer  of  momentum  and  therefore  to 
the  sea-surface  stress.  However,  this  pattern  in  \j/*  may  oscillate  between  \ji*  and  ff* 

throughout  a  period,  and  this  possibility  will  be  explored  later. 

In  Fig.  3.10(a)  and  Fig.  3.10(b),  we  look  at  the  cross-section  perturbation 
temperature  contours  using  the  same  rotated  axes.  In  these  figures,  we  find  that  the 
temperature  exhibits  a  vertical  wavenumber  two  pattern  in  both  cross-sections.  Both  this 
pattern  and  the  streamfunction  patterns  are  similar  to  Schmidt  and  Schumann’s  (1989) 
temperature  and  streamfunction  patterns  produced  by  their  large-eddy  simulation  (LES). 
Their  LES  run  had  an  inversion,  while  our  NDS  does  not  have  an  inversion;  thus  the 
similarity  in  the  cross-sections  is  quite  gratifying. 

To  complement  the  contoured  streamfunction  diagrams,  we  examine  the  domain 
cross-section  velocity  vectors  whose  horizontal  components  are  u*  and  v*.  Here,  the 
perturbation  velocity  u*  is  given  by  the  curl  of  the  streamfunction  ii|*and  the  perturbation 
velocity  v*  is  given  by  the  curl  of  the  streamfunction * .  Figures  3.11(a)  and  3.11(b) 


y*-across 


(a) 


Fig,  3.8.  Domain  cross-section  of  contoured  streamfunctions:  (a)  \j7*  and  (b)  fj*  in  the  cross-mean-wind 
direction  -across.  The  f|‘  contours  represent  U*  approximately  along  the  mean  wind 

direction  and  the  contours  represent  v*  approximately  perpendicular  to  the  mean  wind 
direction.  The  solid  contours  are  positive  values  and  the  dashed  contours  are  negative  values. 
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Fig.  3.9.  Domain  cross-section  of  contoured  streamfunctions:  (a)  \j?’  and  (b)  fj’  in  the  along-mean-wind 
direction  jc* -along.  The  f]*  contours  represent  w*  approximately  along  the  mean  wind 

direction  and  the  \j/*  contours  represent  v*  approximately  perpendicular  to  the  mean  wind 
direction.  The  solid  contours  are  positive  values  and  the  dashed  contours  are  negative  values. 
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present  the  {v*,w*)  and  velocity  vector  fields  in  the  cross-wind  direction  and 

correspond  to  the  streamfunctions  in  Fig.  3.8(a)  and  Fig.  3.8(b),  respectively.  Figures 
3.12(a)  and  3.12(b)  present  the  (v*,w‘)  and  velocity  vectors  in  the  along-wind 

direction  and  correspond  to  the  streamfunctions  in  Fig.  3.9(a)  and  Fig.  3.9(b), 
respectively.  In  these  figures,  we  see  that  areas  of  divergence  and  convergence  are  offset 
in  the  vertical  suggesting  a  tilt  to  the  circulation.  We  also  find  interesting  small-scale 
circulations  near  the  lower  boundary  in  Fig.  3.1 1(b).  These  circulations  are  responsible 
for  some  downward  transport  of  momentum  to  the  sea-surface. 

The  velocity  vector  cross-sections  display  numerous  strong  updrafts  and 
downdrafts.  We  next  examine  the  dimensional  vertical  velocity  profile  for  the  center  of 
an  updraft  in  Fig.  3.13.  This  graph  indicates  that  in  the  updraft,  the  maximum  vertical 
velocity  occurs  at  z*  =  0.35  and  has  a  typical  magnitude.  This  result  shows  that  by 
controlling  the  vertical  velocity  at  the  lower  boundary,  we  force  the  maximum  vertical 
velocity  up  into  the  middle  of  the  boundary  layer  where  it  actually  occurs.  This  result  is 
consistent  with  convective  boundary  layer  observations  (e.g.,  Lenschow  et  al.  1980)  and 
other  model  simulations.  Therefore,  we  have  solved  Zuccarello’s  (1994)  second  problem 
with  the  model. 

The  model-calculated  winds  affect  the  mean  wind  profile.  The  unrotated, 

dimensional,  horizontally-independent  wind  profile  components,  M-modification  given  by 
r|o^  and  v-modification  given  by  ,  are  plotted  in  Fig.  3.14  and  Fig.  3.15.  These 

steady  profiles  are  independent  of  the  horizontal  position  in  the  domain  and  have  a 
sensible  magnitude.  Both  profiles  show  that  the  ELSE  circulation  is  redistributing  the 
shear  in  the  background  wind  (Haack  and  Shirer,  1992). 
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The  temperature  profile  modification  by  the  circulation,  given  by  ,  is  shown  in 

Fig.  3.16.  This  steady  profile  is  again  independent  of  location  in  the  horizontal.  The 
linear  profile  makes  sense,  as  we  have  heating  only  from  below  and  a  thermally  neutral 
atmosphere,  given  by  Ra  =  -\.  This  plot  helps  justify  the  constant  Ra  approach  of  Haack 
and  Shirer  (1992)  and  others.  The  higher  temperature  at  the  lower  boundary  is  consistent 
with  the  upward  heat  flux  prescribed  by  57.  >  0  in  our  boundary  conditions  (Fig.  2.1). 

In  the  next  step  we  determine  whether  the  model  captures  the  observed  kilometer- 
scale  stress  patterns  seen  in  the  SAR  image  of  Sikora  et  al.  (1995).  We  first  examine  the 
complete  model-calculated  velocity  field  at  the  lower  boundary  in  Fig.  3.17.  This  pattern 
is  obtained  from  the  vector  sum  of  the  circulation  velocity  vectors  and  the  wind  profile 
modification  velocity  vectors,  which  are  expressed  relative  to  the  ocean  current.  The 
velocity  field  displays  a  distinct  three-dimensional  cellular  pattern,  as  expected  from  the 
above  results,  with  a  horizontal  spacing  of  1330  m.  The  prominent  orientation  of  the 
lines  of  cells  is  transverse  to  the  background  wind,  and,  from  the  0  value  above,  is  34° 
west  of  north.  The  shape  of  the  domain  is  chosen  to  match  the  dimensions  given  by  the 
preferred  aspect  ratios  for  easy  comparison  with  the  SAR  image. 

By  similarity  theory,  the  horizontal  velocity  at  the  lower  boundary  is  proportional  to 
the  horizontal  velocity  at  the  sea  surface  (e.g.  Stull  1988).  By  the  standard  drag  law 
parameterization,  the  total  horizontal  wind  speed  is  proportional  to  the  square-root  of  the 
sea-surface  stress  magnitude  (e.g.  Stull  1988).  To  find  this  magnitude,  we  therefore  need 
the  vector  sum  of  the  complete  model-calculated  winds  and  the  background  mean  wind, 
which  also  is  expressed  relative  to  the  ocean  current.  The  resulting  stress  pattern  is 
plotted  in  Fig.  3.18,  with  small  values  shaded  more  and  large  values  shaded  less  in  order 
to  match  the  brightness  of  the  stress  pattern  on  the  SAR  image.  Quite  pleasingly,  the 
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Fig.  3.17.  Planview  of  complete  horizontal  dimensionless  velocity  field  at  the  lower  boundary  z*=0. 

The  shape  of  the  domain  is  chosen  to  agree  with  the  ratio  cLyjci^ ,  the  actual  slope  of  the  y  -axis. 


59 


Fig.  3.18.  Planview  of  stress  magnitudes  at  the  lower  boundary  z  =0.  Small  values  are  shaded  more, 
large  values  are  shaded  less,  in  order  to  match  the  brightness  of  the  stress  pattern  on  the  SAR 
image.  The  domain  shape  is  chosen  as  in  Fig.  3.17. 


60 


stress  pattern  on  this  planview  matches  fairly  well  the  orientation  of  35°  west  of  north 
and  the  spacing  of  approximately  1000  m  observed  in  the  stress  pattern  in  the  SAR  image 
published  in  Sikora  et  al.  (1995). 

3.4  Case  Study  Conclusions 

We  are  quite  pleased  with  the  way  the  model  is  now  running.  The  energetics 
analysis  and  the  incorporation  of  the  new  boundary  conditions  enabled  us  to  identify  and 
correct  some  coding  errors,  as  well  as  to  streamline  the  operation  of  the  model.  Quite 
significantly,  the  model  results  are  consistent  with  results  from  previous  studies, 
including  higher-resolution  LES. 

We  have  several  improvements  to  add  to  the  model  in  the  near  future.  The  addition 
of  a  temperature  inversion  will  cap  the  circulations  more  realistically  than  will  use  of  a 
rigid  lid  upper  boundary  condition  (Laufersweiler  and  Shirer  1995).  We  have  a  separate 
routine  ready  to  be  tested,  that  will  calculate  the  mean-wind  profile  parameters  and 
dy  needed  to  approximate  other  known  profiles,  such  as  an  Ekman  profile,  with  profiles 
obeying  our  boundary  conditions.  We  will  use  these  profiles  to  help  us  compare  our 
preferred  aspect  ratio  and  orientation  results  with  those  obtained  from  stress-free  models. 
We  also  plan  to  calculate  the  sub-BLSE  fluxes,  which  may  help  us  to  better  understand 
the  heat  and  momentum  flux  profiles  given  by  our  model.  The  current  results  are  very 
promising  and  show  that  the  model  is  now  ready  to  be  used  in  further  studies  of  how 
ELSE  affect  the  sea-surface  stress  variability  in  different  boundary  layer  configurations. 
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Appendix  A 

MODEL  ENERGETICS 


We  analyzed  the  energetics  of  the  dimensionless  model  equations  to  identify  all 
energy  sources  and  sinks.  This  analysis  was  essential  as  a  diagnostic  tool.  We  compared 
the  results  from  the  original  equation  of  motion  and  the  vorticity  equation  used  in  the 
model,  to  determine  the  conditions  ensuring  that  there  were  no  new  energy  sources  or 
sinks  introduced  by  use  of  the  vorticity  equation. 

We  began  the  analysis  by  deriving  the  kinetic  energy  (^£)  rate  equation  from  the 
equation  of  motion.  Using  the  available  potential  energy  (AE)  to  KE  conversion  term  in 
the  KE  rate  equation  to  provide  the  definition  for  AE,  we  derived  the  AE  rate  equation 
from  the  thermodynamic  equation.  We  derived  another  KE  rate  equation  from  the 
vorticity  equation.  We  concluded  the  analysis  by  comparing  the  energy  rate  equations, 
justifying  the  energy  sources  and  sinks,  and  identifying  the  boundary  conditions  needed  to 
make  both  equations  the  same. 

We  derived  the  first  KE  rate  equation  by  first  multiplying  the  equation  of  motion 
by  the  velocity  vector  v*  and  then  taking  the  volume  integral  of  the  product: 


jw 


Vp  +PV  V  -/  Pkxv  +PT  k 

ot  ^ 


-(ReV*  +v*)«  Vv*  -Rew* 


(A.1) 
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The  horizontal  aspect  ratios  are  a,  =  and  aj  =  »  we  set  Oj  =  =  1 .  In  the  equations 

below,  we  use  the  notation  in  the  last  term  in  (A.2). 

We  rewrite  the  first  term  so  that  the  time  derivative  is  outside  the  volume  integral. 
The  result  is  a  time  derivative  of  KE,  denoted  KE ,  as  follows: 


V  V  V 


The  result  of  calculating  the  vector  dot  product  of  the  right  side  of  the  equation,  using 
integration  by  parts,  introducing  the  boundary  conditions  (2.34)  -  (2.39),  and  integrating 
over  z*  when  possible,  is: 


+  w 


dz 
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Rew  V,.  -r-^dV , 
dz 


(4)  (5)  (6) 


(A.4) 


where  B  is  the  horizontal  boundary  and  we  use  the  definition  of  V  in  (A.2). 

The  six  terms  on  the  right  side  of  (A.4)  represent  sources  and  sinks  of  energy. 
The  first  term  is  a  boundary  term  due  to  pressure  work.  The  second  term  is  traditionally 
the  KE  dissipation  term  and  is  clearly  an  energy  sink.  The  third  and  fifth  terms  are  also 
boundary  terms.  The  three  boundary  terms  (1),  (3),  and  (5)  are  either  energy  sources  or 
sinks  depending  on  the  sign  of  the  velocity  term  within  the  integral  and  the  sign  of  the 
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velocity  Schramm  constant,  in  term  (3)  (for  this  model  >0).  The  fourth  term  is  a 

convective  energy  source  term  describing  the  conversion  of  AE  to  KE.  We  denote  it  as 
C(AE,  KE) .  The  last  term  (6)  is  a  KE  generation  or  Reynolds  stress  term  representing 

the  dynamic  source  of  energy. 

We  used  the  C{AE,KE)  term  in  (A.4),  C{AE,KE)  =  j  PT*w*dV,  to  derive  the 

AE  rate  equation,  denoted  AE .  We  need  the  same  term,  but  with  the  opposite  sign,  in 
the  AE  rate  equation.  We  derived  the  AE  equation  by  multiplying  the  volume  integral  of 
the  thermodynamic  equation  by  an  appropriate  factor  to  obtain  -C(AE,  KE) .  At  first,  we 

determined  this  factor  to  be  -PT*  /  Ra : 


=V  V*  -ReV*  -Vr*  +Raw*  -v*  •VT 


(A.5) 


However,  looking  at  the  first  term,  we  see  that 


a 


=  AE. 


(A.6) 


This  term  is  valid  only  for  the  case  Ra<0,  since  AE  must  be  nonnegative. 

In  our  model  we  want  the  value  of  Ra  to  be  unrestricted  and  so  we  introduce  a 
correction  term  by  adding  and  subtracting  the  vertical  velocity,  -w*  +  w* ,  to  the 
thermodynamic  equation.  We  then  multiply  the  thermodynamic  equation  by  PT*  to 
obtain: 


J([ 


3-»-=V^r  -ReV  •Vr  +Raw  -w  +w 
dt 


-v‘«vr‘]*(pr)}jv.  (A.7) 
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We  now  get  an  AE  definition  that  holds  for  any  value  of  Ra: 


dV  =  AE. 

J 


(A.8) 


After  integrating  by  parts,  applying  the  boundary  conditions  (2.34)  -  (2.39),  and 
integrating  over  z*  when  possible,  we  get  the  AE  equation: 


dr 
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dV  + 


I 


V 


{Ra  +  \)Prw*dV 
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(4)  (5) 


(A.9) 


The  six  terms  on  the  right  side  of  (A.9)  represent  sources  and  sinks  of  AE.  The 
first  term  is  a  boundary  term  that  depends  on  the  sign  of  the  temperature  Schramm 
constant,  Sj .  We  use  Sj>Q  in  this  model,  so  the  term  (1)  is  an  AE  source.  The  second 

term  is  an  AE  dissipation  term  and  is  clearly  a  sink.  The  third  term  comes  from  the 
original  Rayleigh  term  in  the  thermodynamic  equation  and  our  correction  term.  This  term 
is  an  AE  source  for  Ra>-\  and  a  sink  for  Ra<-\.  The  case  of  Ra  =  -\  corresponds  to 
a  thermodynamically  neutral  atmosphere.  In  term  (4),  we  find  the  C(AE,KE)  term  of  the 
KE  equation,  but  with  the  opposite  sign.  It  represents  a  source  or  sink  of  AE,  but  not  a 
net  source  of  total  energy  KE+AE.  The  last  term  is  another  boundary  term,  and, 
depending  on  the  sign  of  the  vertical  velocity  at  z*  =  0 ,  is  either  a  source  or  sink  of  AE. 
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We  completed  the  energetics  analysis  for  the  original  equations  of  the  problem 
based  on  our  original  boundary  conditions.  Now  we  need  to  analyze  the  vorticity 
equation  used  in  the  model.  For  physical  consistency,  we  must  obtain  the  same 
KE  equation  found  from  the  equation  of  motion  (A.3).  To  get  this  form,  we  recognize 
that  the  velocity  perturbation  in  the  model  (Zuccarello,  1994)  is  v*  =  V  x  A* ,  where  the 
vector  streamfunction  A*  is  the  sum  of  the  two  vector  streamfunctions  Xj/  *  and  r|*: 

v*  =  VxA*  =Vx\i/*+Vxri*.  (A.IO) 


We  then  use  a  vector  identity  and  rearrange  its  terms  to  give: 

A*  •(Vx  v‘)  =  V«(v*  X  A‘)  +  v‘  •(Vx  A*). 


We  then  substitute  v*  =  V  x  A*  into  the  last  term  to  obtain: 

A-.(Vxv‘)  =  V.(v*xA‘)  +  v*.v‘. 


(A.11) 


(A.12) 


So,  by  vector  multiplying  the  vorticity  equation  by  the  vector  streamfunction  A,  we  get 
the  dot  product  of  the  velocity  vector  needed  to  form  KE  ,  together  with  another  term. 

We  multiply  the  vorticity  equation  by  A  and  integrate  the  product  over  the  volume 
V: 


PV'  (V  X  v* )  -  /*P(V  X  (k  X  v‘ ))  +  P(V  X  r*k) 
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The  result,  after  vector  multiplication,  integration  by  parts,  application  of  (A.  10),  (A.  12) 
and  the  boundary  conditions  (2.34)  -  (2.39),  and  integration  over  z  when  possible,  is: 
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z  =0 


(12) 

where  the  tensor  is  the  streamfunction  component  \|/*  when  ^  =  1 ,  Ti*  when  k  =  2, 
and  0  when  k  =  3. 

Now,  we  compare  the  two  KE  equations  (A.4)  and  (A.14).  We  do  not  find  term 
(1),  the  pressure  boundary  term  in  (A.4),  directly  in  (A.14).  However,  we  find  the  next 
five  terms  in  both  (A.4)  and  (A.14).  The  terms  are:  the  dissipation  term  (2);  a  boundary 
term  (3);  the  C(AE,  KE)  term  (4);  another  boundary  term  (5);  and  the  KE  generation  term 

(6).  We  can  find  no  direct  comparison  in  (A.4)  with  the  six  remaining  terms  in  (A.14). 
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Therefore,  these  six  terms  must  come  from  the  pressure  boundary  term  in  (A.4),  since 
these  two  equations  must  be  equal  and  the  pressure  boundary  term  is  the  only  term  not 
found  in  (A.  14). 

It  would  be  much  easier  to  study  the  energetics  of  the  model  if  we  could  eliminate 

•  • 

the  unwanted  sources  and  sinks  of  energy.  In  both  KE  equations  and  the  AE  equation, 
the  unwanted  sources  and  sinks  are  boundary  integrals  with  w*[z*  =  0)  and  A,  (z*  =  0)  in 
them.  By  making  w*(z*  =  0)  =  0  and  A,(z*  =  O)  =  0,  as  explained  in  Appendix  B,  we 
make  these  boundary  terms  vanish.  (Note  that  we  assume  dw*  jdz*  is  finite  at  z*  =  0 , 
which  is  a  realistic  condition.)  With  these  boundary  conditions,  the  energy  rate  equations 
(A.4)  =  (A.  14)  and  (A.9)  become: 


(A.  15) 
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(1)  (2)  (3)  (4) 

We  have  rectified  the  model  energetics.  In  the  KE  equation  the  terms  are:  KE 
dissipation  (1),  a  boundary  KE  sink  term  (2)  (since  >  0),  the  C(AE,KE^  term  (3),  and 

the  KE  generation  term  (4).  Similar  terms  are  in  the  AE  equation:  a  boundary  AE  source 
term  (1)  (since  5^  >  0 ),  the  AE  dissipation  term  (2),  the  AE  generation  term  (3),  and  the 

C(AE,KE)  term  (4).  Note  that  for  Ra  =  -\,  the  boundary  layer  flow  is  driven  only  by  heat 
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flux  from  the  bottom  and  modified  only  by  shear.  We  examine  this  special  case  in  this 
thesis. 


72 


Appendix  B 


MODMCATION  OF  THE  VERTICAL  BASIS  FUNCTION 


The  energetics  analysis,  discussed  in  Appendix  A,  resulted  in  boundary  conditions 
A‘(z*  =  O)  =  0 ,  which  imply  w*(z‘  =  O)  =  0,  that  supplement  the  original  Monin- 

Obukhov  boundary  conditions  (2.34)  -  (2.39).  After  discovering  the  need  for  these  new 
vertical  boundary  conditions,  we  did  not  want  to  completely  rederive  the  model. 
Therefore,  we  reformulated  the  vertical  velocity  basis  functions  to  satisfy  both  the  new 
and  the  old  boundary  conditions,  and  then  we  changed  the  program  to  account  for  these 
changes.  This  appendix  explains  the  vertical  basis  function  reformulation  and  the 
modifications  to  the  existing  program. 

B.l.  Vertical  Basis  Function  Reformulation 

We  start  by  looking  at  the  condition  A*(0)  =  0,  and  then  extending  it  to  the 
condition  w*(0)  =  0.  We  represent  the  velocity  vector  v*  as  the  curl  of  a  vector  stream- 
function  A* ,  which  is  the  sum  of  two  vector  streamfunctions  *  and  r|* .  Therefore,  we 
may  express  the  velocity  perturbation  v*  as 

V*  =  VxA*  =  Vxv|/*-i-VxTi*.  (B.l) 

We  use  a  truncated  Fourier  expansion  to  represent  the  dependent  vector  streamfunctions 
\j/*  and  ri*.  Accordingly,  the  vector  streamfunction  A*  is 
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4  4 


p=0  q=\ 


■  (B.2) 


+ 


where  the  only  dependence  on  z*  is  in  the  vertical  basis  function,  h^[z')  ■  If  we  want 
A*(z*  =  O)  =  0 ,  then  we  must  force  \{0)  =  0 ,  for  each  of  the  vertical  wavenumbers,  q. 

The  vertical  basis  function  is 


/j^(z*)  =  cosG3^Z*-^„05^sinG5^Z*  for  ^  =  1,2,3,4  , 


(B.3) 


where  G5 , ,  03  2 ,  G5  3 , 05  4  are  the  first  four  positive  solutions  of 

5„,03^=cot03,.  (B.4) 


The  vertical  basis  function  values  at  the  vertical  boundaries  are: 


/J,(1)  =  0  (B.5) 

/t4(0)  =  l.  (B.6) 


Using  (B.6),  we  find  a  very  simple  way  to  make  the  new  vertical  basis  functions  vanish 
that  leads  to  A*(z*  =  O)  =  0 .  We  define  new  vertical  basis  functions  ,  for 

=  1, 2, 3 ,  by  taking  appropriate  differences  of  the  old  basis  functions  : 


Kew\ 

Kewl 

Kew3 


lA 

“  Km  ~  Kidi 

(B.7) 

“  Kidi  ””  Kid'i 

(B.8) 

^  Km  ~  KidA  • 

■  (B.9) 

The  three  new  vertical  basis  functions  reduce  the  number  of  vertical  wavenumbers  in  the 
vorticity  representation  from  four  to  three.  We  use  (B.3)  to  define  the  old  basis  functions 
and  then  use  these  to  obtain  the  new  basis  functions  (B.7)  -  (B.9).  We  use  these  new 
functions  to  plot  the  vertical  velocity  profiles,  as  shown  in  Fig.  B.l. 

Once  put  into  a  new  Galerkin  expansion,  the  resulting  new  function  satisfies  the 
boundary  condition  A*(z*  =  O)  =  0 : 


(x* ,  y*  ,z*  ,r*)  =  (r*  (x* ,  y*  (z*  )}i 

P^0q=\ 

+{fi  pg'  ‘  y^Sp  y)hne^^'  {z  )}j  +  Ok]  . 


(B.IO) 


Substituting  (B.7)-(B.9)  in  (B.IO)  yields 


A  *  (x* ,  y*  ,z* ,  t  * )  =  X  X  [{'K  (O^^^Sp  {x  ,  y*)(Ku^'  (z* )  -  (z*  ))}i 

p^Oq'^l 

+{^  P,'  (^*  )(Kw  (^* )  -  Kw.x  (^* ))} j  +  Ok] . 


(B.ll) 


At  z*  =  0 ,  the  old  basis  functions  =  1  by  (B.6)  and  therefore  their  differences  in 
(B.ll)  cancel.  Summing  over p  and  q'  results  in  A*(z*  =  O)  =  0 . 

When  we  carry  out  the  curl  of  A*  in  (B.l  1),  we  get  the  new  velocity  component 
equations.  The  vertical  velocity  equation  is: 


VELOCITY  BASIS  FUNCTION  (DIMENSIONLESS) 


Fig.  B.l.  Velocity  vertical  basis  function  (vertical  velocity)  profiles  (B.7)  -  (B.9)  for  the  three  new 
wavenumbers. 


Summing  over  p  and  q' ,  we  get  w*(0)  =  0 ,  which  is  our  other  new  boundary  condition. 
This  accomplishes  our  first  goal  of  making  w*(0)  =  0  and  A‘(0)  =  0 . 

B.2.  Modification  of  Existing  Program  to  Incorporate  New  Basis  Functions 

To  account  for  the  new  boundary  conditions,  we  adjust  the  existing  model  by 
rewriting  the  portions  of  the  code  involving  the  vertical  velocity  basis  function  .  When 

we  define  the  new  vertical  basis  function  ,  we  reduce  the  streamfunction  vertical 

wavenumbers  from  four  to  three.  This  result  reduces  the  number  of  model  equations  and 
coefficients  from  60  to  50.  We  now  integrate  the  model  with  a  50-element  amplitude 
coefficient  vector  y.  Therefore,  as  for  the  60  coefficient  models,  we  construct  each  of  the 
matrices:  the  temporal  derivative  matrix  A,  the  linear  matrix  B,  and  the  y-dependent 
nonlinear  matrix  C,  as  50  x  50  matrices. 

Ay  =  By  +  C(y)y,  (B.13) 

In  the  linear  portion  of  the  program,  we  introduce  the  new  vertical  basis  function 
definitions  (B.7)  -  (B.9)  into  the  appropriate  integrals.  This  process  is  straightforward 
and  produces  the  50  x  50  matrices  A  and  B.  It  is  a  little  more  complicated  for  the  y- 
dependent  nonlinear  matrix  C.  Zuccarello  (1994)  derives  the  nonlinear  portion  of  the 
model  using  DERIVE  to  find  the  integrals  of  the  products  of  the  vertical  basis  functions. 
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He  defines  these  integrals  using  two-  and  three-dimensional  arrays  in  the  nonlinear 
subroutine.  To  avoid  rederiving  with  (B.7)  -  (B.9),  the  equations  used  to  define  the 

arrays  for  C,  we  found  we  could  index  and  appropriately  difference  these  arrays  to  get  the 
new  integral,  that  are  composed  of  the  products  of  the  new  vertical  basis  functions  . 

The  result  is  a  50  x  50  matrix  C.  Thus  the  model  is  strictly  a  50  coefficient  system,  and  is 
integrated  numerically  as  such. 


